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The Ewald method was originally invented to compute the Madelung 
constant. In this paper we consider a lattice whose sites are associated with 
an arbitrary potential junction. The "charge," or the scale factor for these 
potential functions, need not be the same at each site. We consider the evalua- 
tion of the resulting lattice sum at an arbitrary point, not necessarily at a 
lattice site. The method involves two generalizations over previous work: 
(1 ) the displacement of the origin off a lattice site and {2) the handling of 
arbitrary periodic charge distributions by decomposing such distributions 
into simpler ones involving only -\-q and —q. The method should prove 
particularly useful for evaluating the expansion coefficients of the crystalline 
potential when this potential is expanded in the usual spherical harmonic 
series. 

The problem of summing slowly converging series is an old one. One 
physical context in which the problem has been widely studied is the 
calculation of the potential due to an ionic crystal lattice. The methods 
of Madelung 1 and Evjen 2 depend on collecting ions into neutral groups. 
The convergence obtained in this way, however, is conditional: that is, 
the result depends on the way in which the neutral groups are chosen. 
Ewald 's 3 method, which hinges on doing part of the summation in 
reciprocal space, gives rapid convergence and the limit is unique. 
Subsequent discussions 4 " 1 ' of this topic have been extensions and general- 
izations of these methods. This work too is an extension of the Ewald 
technique. In particular it is a generalization of the approach taken by 
Nijboer and DeWette. 9 

For purposes of orientation, we summarize the basic philosophy of 
the Ewald method. Suppose we have a function tp{r) such that the series 

S=Z<p(T„) (1) 

n-1 

427 
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is slowly converging. The symbol 2^» is to be understood as a shorthand 
for 2Z»i 2]» 2 2Z«3 • ^ represents independent summation on all three 
components of the vector r. We now construct some function ^(r), 
which falls off rapidly with r, and its partner 

/(r) = 1 - g(r), (2) 

which rapidly approaches unity as r increases. We now write 

S = Z <p(rn)Q(r n ) + E *»(r»)/(r„). (3) 

The first sum converges rapidly, because of g. The second sum con- 
verges like <p, i.e., slowly. Its Fourier transform, however, will converge 
rapidly. In fact the more slowly this sum converges the more rapidly 
will its transform converge. To complete the argument we need Parse- 
val's theorem : 
If 

$(h) = f exp (i2irh-rV(r) dr (4a) 

and 

/.'(h) = f exp (i27rh-r)/(r) dr (4b) 

then 

f $(h)F*(h) dh = f ^(r)/*(r) dr (4c) 

where the symbol * denotes "complex conjugate." The formal passage 
from sums to integrals can be accomplished by means of Dirac delta 
functions 5(r — r„), as we shall see below. Thus ParsevaFs theorem 
guarantees that the summation in transform space yields the same 
result as the summation in the original coordinate space. 

We now apply this scheme to the calculation of the potential due 
to an ionic lattice. To begin with, we consider what we shall call a 
"primitive" lattice. Such a lattice is generated from primitive transla- 
tions Ci , c 2 , c 3 in such fashion that 

3 

r„ = 2Z nfii , (5) 

with n\ , n 2 , n 3 taking on independently all integer values from -» to 
oo ; and in addition there is associated with each lattice point r n a charge 

q n = q (-l) ni+n2+n3 (6) 
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where </o is some constant. A typical primitive lattice is NaCl, in con- 
trast, for instance, to CaF 2 , which does not obey (C). 
We define reciprocal vectors h, by the usual relation 

h.-C; = Sij (7) 

and the reciprocal lattice as the aggregate of points 

h„ = 5><hi (8) 

where the n's again run from — <» to a> . It is trivial to show that if 
these two lattices are represented respectively as J^„5(r — r n ) and 
2^„5(h — h„), then the reciprocal lattice is simply the Fourier trans- 
form of the coordinate lattice, the Fourier transform being understood 
as in (4a) and (4b). If we define the special reciprocal lattice vector 

k = i(hi + h, + h,), (9) 

then (6) can be rewritten 

q„ = q exp (i2rk-r n ). (10) 

In addition to q„ , we associate with each lattice point a function 
ip(r). For the present we place no restriction on <p{i), except that it 
possess a Fourier transform. Of course there would be no practical 
motivation for the calculation unless <p(r) fell off slowly with r. We 
wish to sum the contribution of all the p's at some arbitrary point R: 

S = E v(r n - R) exp (i2wk-T n ). (11) 

n 

To change the sum into an integral, as required for the eventual 
application of (4c), we define 

w(t) = exp (*2rk-r) £ 5(r - r„), (12) 

n 

so that 

S = J w(t)tp(t - R)dr. (13) 

In exact analogy to (2 ) and (3 ) we can break £ into two parts : 

8 = f w(r)<p(r - R)</(r- R)dr 

(14) 
+ f u>(rV(r - R)/(r - R)rfr. 

The first integral in (14) corresponds to the first sum in (3) : 
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2>(r„ - R)g(r„ - R)(-l)" 1+ " 2+ " 3 . (15) 

n 

The second integral we wish to evaluate in the conjugate domain. For 
brevity we define 

f(r) - *>(r)/(r) (16) 

tf(h) = [ exp (*2irK-r) *(r) dr. (17) 
Then the Fourier transform of <p(r — R)/(r - R) is given by 

J exp (i27rh-r) r£(r - R) dr = exp &rii»K ¥(h). (18) 

•'—00 

The transform of w(r) is easily evaluated (see Ref. 9) and is given by 
f exp (fiirh-r) to(r) dr = i £ «(h + k - h„), (19) 

J— oo V c n 

where v c equals Ci-c 2 X c 3 , or the volume of the coordinate unit cell. 
By Parseval's theorem [(4a), (4b), (4c)], the second integral of (14) 
now becomes 



- 2 exp \i2ir{h n - k) -R] *(h n - k). 



(20) 



This completes the essential derivation, since 8 is now expressed in 
terms of the two sums (15) and (20), both of which converge rapidly. 

We consider some special cases. If R = — that is, if we are sitting 
at a lattice point — we presumably will want to exclude the contribu- 
tion of that point itself. If <p(0) is finite, the contribution can be sub- 
tracted outright. If <p(0) diverges, as is usually the case, one must be 
clever about picking the functions g and / so that ^ (0) = f(0)<p (0) does 
not diverge. One then simply omits from the sum (15) the term for 
n = 0, and subtracts from the sum (20) the quantity ^(0). Note that 
one subtracts ^(0), not \K0). The Ewald calculation is obtained in 
this way, if one takes 



*(r) = ur 

g(T) =Erfc (|r|) 



(21a) 

(21b) 

/(r) = Erf (|r|). (21c) 

We note that if <p is real (for example, any central potential ) and if 
we pick g real, then ^ will be real also. The function w is both real and 
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symmetric because of our particular definition of the reciprocal vector 
k. It follows that the sum S as defined in (14) is real. But if we look 
at the partial sum (20), this reality is not at first sight apparent, since 
the R can be chosen arbitrarily. But \p real implies that ^ is Hermitian: 
*(-h n 4- k) = ¥*(h„ - k), and clearly exp [z'27r(h„ - k)-R] is Her- 
mitian. Again, because of the peculiar choice of k, the arguments 
± (h„ — k) are bound to occur in pairs. Since the sum of a function 
and its complex conjugate is real, the sum (20) is always real, which 
we may emphasize by rewriting it: 

9 + 

- Re 2 exp U"27r(h„ - k) -R] *(h„ - k). (22) 

V e n 

+ 

Here £ means that we sum only over half the space. This can be 
accomplished by summing n 3 , for instance, from 1 to °o instead of 
from — oo to co . 

If ^(r) is symmetric, ^(h) will be real, and in the sum (22) we will 
obtain only cosine terms. Conversely, if ^(r) is antisymmetric, the re- 
ciprocal sum will contain only sine terms. Again this is independent 
of the choice of R. 

We now also see clearly why the Ewald method "works." The con- 
vergence difficulties with the series S of (1) concern its asymptotic 
behavior. But this behavior is related to the behavior at the origin of 
the series in reciprocal space. By means of the vector k, we guarantee 
avoidance of the origin in reciprocal space, regardless of any other 
conditions in the problem. 

The sums that most frequently occur in practice are related to the 
expansion of the crystalline potential in spherical harmonics: 

V(x) = E £ C, m | x \ l F,._ m (ft, , Vs ) (23) 

C "« = PI ? ^ I F " I"' _1 F| - (d <» ' *J (24) 

Yi, m (6#) = —. . : U- f e P lm {6). (25) 

L 4tt (I + \m |)! J 

The notation is well known and conventional. Our definition of the 
spherical harmonics Yi m implies Y [m * = F^_ m . Also Yi m (t — 6, -k + <p) = 
( — l)'Yi m (0,tp). The evaluation of the crystal sums Ci m has been 
discussed by Nijboer and DeWette. 9 In our notation, <p(r) here cor- 
responds to r~ l ~ l Yi m (6,<p). Nijboer and DeWette's choice of g is the 
incomplete gamma function 11 
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(J ( T ) = r(» + l,irr 2 )/Y(n + I) (26a) 

f(t) = 1 - g(r) = y(n + l,*r*)/T(n + I). (2Gb) 

They solve the problem subject to the restrictions that (a) the po- 
tential is expanded about a lattice point, i.e., R = 0, and (b) the lattice 
is primitive, in the sense of (5), (6) and (10). Our discussion has made 
it obvious how to remove the first restriction. Generalizing their re- 
sult, via (20), 

c, m (r-R) = 2 7TTr(r+!) 

•[£ I U - R \~**TQ + J,» | r n - R | 2 )F, m (0r n -R , ^„-r) 

n 

•exp (Oik-*) (27) 

+ <V"V*E exp (i27r(h n - k)-R)| h„ - k l'" 2 

n 

• r(l,7r I h n - k | 2 ) X Y lm (d bn . k , vh„-k)]. 

The next generalization of the procedure lies in its application to 
nonprimitive lattices. Such application will clearly be possible if an 
arbitrary lattice can be decomposed into a sum of component primitive 
lattices. We illustrate what we mean by a one-dimensional example. 
Consider a one-dimensional lattice with charges distributed as follows : 

L = 2 1 -3 2 1 -3 2 1 -3 • ■■ . 

Thus we have qi = 2, q 2 = 1, q 3 = 0, q 4 = -3. If the distance between 
successive g's is 1 distance unit, then the basic periodicity is 4. We note 
that £? = °> since we must have a neutral lattice, and that zero is 
itself an allowable q value. Now consider the following sequences of 
numbers of periodicity 4: 

T — 3, _i i — i . . • 
iJ l — 2 2 2 2 

We can represent L as the sum L = 2L X + \/2 L 2 + 2 \/2 L 3 . We 
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note that Lj , L> , Li are all primitive since they fulfill both (5) and the 
two equivalent equations ((i) and (10). In terms of (5), for L X} | c | = ] ; 
for Lo and L 3 , | c | = 2, but their origin of coordinates is shifted by 
one unit. 

We can define a dot product for the L's in the following sense. Sup- 
pose we have L A = ]Ct</« U) ^( r — *.) and L B = S<ff< w 3(r — r,). 
Then L A -L a = ]£«ff< ff< ( *\ where the sum runs over a unit cell. 
We note that, in our example, the components L< are orthonormal, 
and that the projection of L on each L, can be found by taking the dot 
product. 

We also observe that we have three components, which is just enough 
to account for 4 numbers which are subject to the one constraint 
^T, Q> = 0- Out three components L, "span the space" of L, and it is 
clear in general that if L consists of a periodic sequence of P numbers 
whose sum is zero, we shall need (P — 1 ) primitive components. 

The question is: Is it possible to decompose an arbitrary periodic 
sequence L into orthonormal primitive components? Can one devise 
an algorithm for finding these components? Can one do this in three 
dimensions? In considering these questions, we have collaborated with 
Dr. R. L. Graham, and we are particularly indebted to him for pointing 
out the great simplification that results if one confines oneself to 
sequences of periodicity 2" in each dimension. 

Consider a linear sequence L of periodicity 2". If c is the primitive 
translation of L, we define a = c/2". The basic vectors for generating 
primitive component lattices are b ( " = a, b <2) = 2b (1) , ••• b (n) = 2b (n-1) . 
Each b <,} for i > 1, will generate a set of primitive lattices differing 
only in their choice of origin. There will evidently be n such sets of 
components. Primitive component lattices containing 2'" nonzero entries 
per unit cell will have a basis vector of length | c |/2" -f " and will have 
2 n ~ m possible shifts of origin. We note that £ m =r 2 n ~" = 2" - 1, which 
is the correct total number of components. The number of primitive 
components of the same periodicity but with shifted origins doubles 
every time the length of the generating basic vector b, doubles, and of 
course the number of nonzero entries per unit cell halves at the same 
time. The set of origin positions {R,j associated with b, is clearly the 
set of all translations R such that R; — R,- ^ b,- . This includes the set 
{R,_i} plus a new set formed by adding b,_i to all R in {R,-_i}. 

All the primitive component lattices are orthogonal to each other, 
in the sense that L t Lj = in . Within each ''phase-shifted" set, each 



434 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1965 

component has numbers in locations where all the other components 
have zero, so that the members of such a set are obviously orthogonal. 
Now consider two sequences L, and L<_i belonging to sets generated 
by b, and b»_i . Either L,_i will have numbers only where L, has zero. 
Otherwise, each alternate number in L,_i will have zero as a partner 
in Li ; the remaining numbers in L,_i all have the same sign, but their 
nonzero partners in L, have alternating signs. Hence once again the 
dot product is zero. The same argument applies clearly not only to 
members of contiguous sets L, and L«_i , but to members of any two 
sets, Li and Lj , i ^ j. To produce not merely orthogonality, but ortho- 
normality over the unit cell, the normalization factor, or q in (6) and 
(10), must clearly be 2~ m/2 for a component with 2 m nonzero entries. 

In the example we have given, n = 2, and all the above arguments 
can be seen to be rather trivially verified. 

Extension of the preceding to two dimensions is straightforward. 
We consider a two-dimensional array, doubly periodic with periodicity 
2" X 2". The primitive translations Ci and c 2 carry any q into the cor- 
responding q in another cell. Within each cell, the different q's are sepa- 
rated by multiples of Ci/2" and c 2 /2 n , and we shall call these vectors 
a! and a 2 . As before, we define a set of components of equal periodicity 
by defining the basic vectors which will generate the primitive lattice. 
The translations by which components within a set differ we denote 
by R. The algorithm for producing a complete orthononnal set of 
components is simple: 

bi <n = &i (28a) 

b 2 (1) = a, (28b) 

b « = a, + a 2 (28c) 

b 2 (2) = a, - a s (28d) 

b, , " +2) = 2b/' (28e) 

b, ( " +2) = 2b (,,) (28f) 

g ci) = 2- (29a) 

ff W = g'- 1 ' V2 (29b) 

jR u) } = (30a) 

( R (»+»j = | R oo» + { R <»> + bl (»)j # (30b) 
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The set labels appear as superscripts in (28-30). The vectors b 
are the basis of a primitive lattice, and {R} indicates the set of all 
translations R which relate primitive lattices having the same basis. 
Equation (30b) says that the translations for set n include all the 
translations for set n — 1, plus all translations formed by adding any 
one of the vectors b" to all the R's of set n — 1. Again all components 
are orthonormal over the unit cell, so that the decomposition of any 
2" X 2" dimensional cell can simply be found by taking dot products. 
Normalization is insured by our definition of g n , and orthogonality 
becomes clear from the same type of reasoning as in the linear case, 
which we shall not repeat here. 

In three dimensions, we consider an array triply periodic with perio- 
dicity 2" X 2" X 2". We define Ci , c 2 , c 3 , and ai , a 2 , a 3 in analogy with 
the two-dimensional case. We present the basic vectors for the primitive 
lattices, plus the associated translations: 



w° 


= a 3 


(31a) 


W l) 


= a 2 


(31b) 


b:," 1 


= a 2 + ai 


(31c) 


b,« 


= a 3 + a 2 


Old) 


b 2 (2) 


= -a 3 + a 2 


(31e) 


V 2) 


= ai 


(31f) 


bi w 


= ai + a 2 + a 3 


(31g) 


b 2 (3) 


= a, - a 2 + a 3 


(31h) 


b 3 (3) 


= ai + a 2 - a 3 


(31i) 


bi (*+3) 


= 2b 1 (n) 


(31j) 


b 2 ( " +3) 


= 2b 2 ( '° 


(31k) 


b 3 ( " +3) 


= 2b 3 (n) 


(311) 


<z (1) 


_ f)-3 n/2 


(32a) 


ff C«-K> 


= q"V2 


(32b) 


(R (1) ) 


= 


(33a) 


| R (.+i>j 


= {R n } + {R (n) + W B) }. 


(33b) 
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The choice of b's is not unique, and we have given a set that seems to 
bear a maximum analogy to the two-dimensional case. 

Many lattices can be represented as 2" X 2" X 2" dimensional ar- 
rays. We note that the primitive translations Ci , c 2 , c 5 need not be 
orthogonal. Even if a lattice does not lend itself to such a representation, 
by choosing a fine enough grid (with a correspondingly large number of 
zero charges), one can approximate any lattice to any desired degree of 
accuracy.* 

We now return to the problem of the lattice sum. We will obtain a 
sum of the type (20), or more specifically of the type (27a), for each 
primitive component. For a particular set of components which are 
translated from each other by { R„} , it clearly makes no difference whether 
we sum the contribution of each component at the origin or the contri- 
bution of any one component at the points given by {R n } . (The difference 
between shifting the function and shifting the coordinate system is 
merely a conceptual one.) Hence the shift vectors R„ correspond to 
the position vectors R in expressions (20) and (27a). 

In the present approach, we have completely split the geometrical 
character of the lattice from the charges assigned to each lattice point. 
This suggests the possibility of computing the lattice sums arising from 
all the primitive components of commonly occurring grids, once and 
for all. The problem of computing a particular lattice sum would then 
involve only the decomposition of the given lattice into primitive com- 
ponent lattices, whose contribution to the sum would already be known. 
Work along this line is in progress. 
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